home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / coulomb_bound.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  3.7 KB  |  120 lines
  1. /* specfunc/coulomb_bound.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_exp.h>
  26. #include <gsl/gsl_sf_gamma.h>
  27. #include <gsl/gsl_sf_pow_int.h>
  28. #include <gsl/gsl_sf_laguerre.h>
  29. #include <gsl/gsl_sf_coulomb.h>
  30.  
  31. #include "error.h"
  32. #include "check.h"
  33.  
  34. /* normalization for hydrogenic wave functions */
  35. static
  36. int
  37. R_norm(const int n, const int l, const double Z, gsl_sf_result * result)
  38. {
  39.   double A   = 2.0*Z/n;
  40.   double pre = sqrt(A*A*A /(2.0*n));
  41.   gsl_sf_result ln_a, ln_b;
  42.   gsl_sf_result ex;
  43.   int stat_a = gsl_sf_lnfact_e(n+l, &ln_a);
  44.   int stat_b = gsl_sf_lnfact_e(n-l-1, &ln_b);
  45.   double diff_val = 0.5*(ln_b.val - ln_a.val);
  46.   double diff_err = 0.5*(ln_b.err + ln_a.err) + GSL_DBL_EPSILON * fabs(diff_val);
  47.   int stat_e = gsl_sf_exp_err_e(diff_val, diff_err, &ex);
  48.   result->val  = pre * ex.val;
  49.   result->err  = pre * ex.err;
  50.   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  51.   return GSL_ERROR_SELECT_3(stat_e, stat_a, stat_b);
  52. }
  53.  
  54.  
  55. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  56.  
  57. int
  58. gsl_sf_hydrogenicR_1_e(const double Z, const double r, gsl_sf_result * result)
  59. {
  60.   if(Z > 0.0 && r >= 0.0) {
  61.     double A = 2.0*Z;
  62.     double norm = A*sqrt(Z);
  63.     double ea = exp(-Z*r);
  64.     result->val = norm*ea;
  65.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) * fabs(Z*r);
  66.     CHECK_UNDERFLOW(result);
  67.     return GSL_SUCCESS;
  68.   }
  69.   else {
  70.     DOMAIN_ERROR(result);
  71.   }
  72. }
  73.  
  74.  
  75. int
  76. gsl_sf_hydrogenicR_e(const int n, const int l,
  77.                         const double Z, const double r,
  78.                         gsl_sf_result * result)
  79. {
  80.   if(n < 1 || l > n-1 || Z <= 0.0 || r < 0.0) {
  81.     DOMAIN_ERROR(result);
  82.   }
  83.   else {
  84.     double A = 2.0*Z/n;
  85.     gsl_sf_result norm;
  86.     int stat_norm = R_norm(n, l, Z, &norm);
  87.     double rho = A*r;
  88.     double ea = exp(-0.5*rho);
  89.     double pp = gsl_sf_pow_int(rho, l);
  90.     gsl_sf_result lag;
  91.     int stat_lag = gsl_sf_laguerre_n_e(n-l-1, 2*l+1, rho, &lag);
  92.     double W_val = norm.val * ea * pp;
  93.     double W_err = norm.err * ea * pp;
  94.     W_err += norm.val * ((0.5*rho + 1.0) * GSL_DBL_EPSILON) * ea * pp;
  95.     W_err += norm.val * ea * ((l+1.0) * GSL_DBL_EPSILON) * pp;
  96.     result->val  = W_val * lag.val;
  97.     result->err  = W_val * lag.err + W_err * fabs(lag.val);
  98.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  99.     if (stat_lag == GSL_SUCCESS && stat_norm == GSL_SUCCESS) {
  100.       CHECK_UNDERFLOW(result);
  101.     };
  102.     return GSL_ERROR_SELECT_2(stat_lag, stat_norm);
  103.   }
  104. }
  105.  
  106. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  107.  
  108. #include "eval.h"
  109.  
  110. double gsl_sf_hydrogenicR_1(const double Z, const double r)
  111. {
  112.   EVAL_RESULT(gsl_sf_hydrogenicR_1_e(Z, r, &result));
  113. }
  114.  
  115.  
  116. double gsl_sf_hydrogenicR(const int n, const int l, const double Z, const double r)
  117. {
  118.   EVAL_RESULT(gsl_sf_hydrogenicR_e(n, l, Z, r, &result));
  119. }
  120.